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is  used  in  the  area  where  the  flow  is  highly  nonlinear  but  still  subcritical. 
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already  been  published»j[n  the  scientific  literature  (references  1 and  2). 
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FLOW  FIELD  MATCHING 


I.  INTRODUCTION 


In  transonic  flow  problems  with  subsonic  freestream  and  infinite  flow  field, 
discretized  solution  methods  which  by  nature  solve  the  equations  only  over  a finite 
domain  are  not  in  themselves  adequate.  The  influence  on  the  near  field  solution  of 
the  "far  field",  which  extends  to  infinity,  must  be  accounted  for.  The  present 
work  has  explored  a matching  method  which  has  proved  to  be  superior  to  the 
asymptotic  expansion  commonly  employed  to  accomplish  this. 

In  the  problem  under  consideration,  the  equations  reduce  to  a linear,  elliptic 
Prandtl-Glauert  equation  in  the  limit  at  infinity.  That  equation  is  soluble  on 
infinite  regions  by  methods  of  superposition  using  "panels"  of  Green's  function 
integrals.  The  simplest  matching  method  introduces  a matching  boundary  in  the 
flow  field  well  outside  the  supersonic  region.  Outside  the  boundary,  a Prandtl- 
Glauert  solution  is  found,  while  inside  a nonlinear  transonic  solution  is  found 
such  that  both  value  and  nonnal  derivative  of  the  velocity  potential  match  across 
the  boundary. 

Two  more  refined  methods  were  tried.  In  one,  a Poisson  term  was  added  to  the 
Prandtl-Glauert  equation  to  make  it  a second-order  approximation  to  the  nonlinear 
equation.  In  the  other,  a "middle  field"  of  large  finite  elements  was  introduced 
between  the  panels  and  the  (finite  difference)  inner  field,  with  matching  conditions 
as  before. 

In  order  to  avoid  duplication  of  effort,  the  present  report  details  the 
finite  element  results  while  only  outlining  results  that  have  already  been  published 
in  our  papers  and  previous  quarterly  reports. 

II.  THEORY 


The  transonic  small  disturbance  equation  (TSDE)  used  in  the  finite  difference 
and  finite  element  areas  is  in  3D, 


wkert  • 


J 
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The  2D  equation  is  the  same  except  the  y-  derivative  term  is  missing.  The  first 
order  (FO)  far  field  equation  is  the  Prandtl-Glauert  equation 
(2)  • 

The  second  order  (SO)  Poisson  correction  term  used  in  the  2D  second  order  method 
satisfies  . 

(3)  \ ^ ^ X ^ 

where  is  a first  order  solution. 


The  finite  difference  method  used  was  conservative  over-relaxation  of  the 
Murman  type.  The  outer  flow  is  solved  by  source  panels  (see  ref.  4)  along  the 
matching  boundary,  the  solution  being  updated  every  10  to  80  sweeps.  In  the  second- 
order  2D  method,  the  equation  being  solved  for  is  linearized  about  a new  value 
of  about  once  every  five  first  order  updates. 

Two  methods  of  matching  were  tried:  Revising  the  outer  flow  to  partly  eliminate 
the  normal  derivative  discontinuity  between  inner  and  outer  regions  (an  "underelaxation" 
factor  of  .5  being  theoretically  optimal  here),  and  adding  sufficient  source  strength 
to  cancel  this  discontinuity  if  it  were  superimposed  on  both  inner  and  outer  flows. 

The  latter  method  is  simpler  and  easier  to  justify  (ref,  1),  and  both  seem  to  be 
about  equivalent  in  practice. 

For  more  details  on  these  methods  see  refs.  1,  2,  3.1,  3.2,  3.3.  and  Figures  1 and  2. 

The  extension  of  the  second-order  method  to  3D  proved  to  involve  really  unwieldy 
influence-coefficient  calculations,  and  its  stability  was  called  somewhat  into  question 
by  the  later  2D  results  (ref.  3.3).  Therefore  we  decided  (ref.  3.4)  to  introduce  a 
middle  field  with  large  finite  elements  satisfying  an  approximation  to  (1).  The  finite 
element  (FE)  region  is  highly  nonlinear  but  still  subcritical,  and  the  usual  elliptic 
variational  approaches  are  therefore  usable  (see  ref.  5). 

Figure  3 shows  a schematic  example  of  a 3-D  midfield  with  finite  elements.  Both 
potential  and  normal  derivative  continuity  are  enforced  on  the  interface  boundary  between 
FD  and  FE  regions  by  overlapping  by  one  mesh  width.  Since  the  same  equation  is  solved 
in  both  regions,  the  overlapping  gives  enough  resolution  for  the  continuity  requirements. 
Dirichlet  conditions  are  used  on  a11  sides  of  the  FE  region  boundary,  and  the  overlapping 
updates  every  cycle  result  in  convergence  of  the  normal  derivative  jump  to  zero. 

Figure  4 diagrams  the  ordering  by  which  the  three  flow  fields  are  brought  to 
simultaneous  convergence  (see  ref.  3.5). 
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For  the  nonlinear  TSDE,  the  finite  element  formulation  introduces  a variational 
functional  which  will  be  minimized  by  the  first  variational  principle.  Only 
rectangular  elements  are  used,  which  allows  the  potential  to  be  approximated  simply 
by  using  separation  of  variables.  The  finite  element  formulation  is  reviewed  in 
Appendix  A. 

A finite  element  code  was  developed  for  both  trilinear  and  triquadratic  elements. 

The  basic  shape  functions  are  defined  and  compared  in  Appendix  B.  The  trilinear 

function  is  continuous  (C°)  while  the  triquadratic  spline  is  slope  continuous  (C^). 
However,  the  resulting  matrix  for  the  quadratic  elements  has  width  five  along  the 

diagonal,  while  it  is  tridiagonal  for  the  linear  elements. 


III.  RESULTS 

As  described  in  Ref.  2 and  3.1,  the  first-order  and  second-order  methods  were 
very  successful  in  a simple  supercritical  non-lifting  case.  Table  1 shows  a cost 
and  quality  comparison.  Those  results  marked  by  asterisks  were  judged  to  show  no 
significant  errors  due  to  the  far  field  approximation. 


Table  1.  Computer  Times  and  Iterations 


Box 

Far-Field 

CP  Sec 

Sweeps 

FO  It 

SO  It 

Large 

Freestream 

*266 

612 

/ 

/ 

Medium 

Klunker 

114 

507 

/ 

/ 

Medium 

FO 

*136 

597 

10 

/ 

Medium 

SO 

*158 

682 

12 

3 

Small 

Klunker 

34 

238 

/ 

/ 

Small 

FO 

82 

542 

10 

/ 

Small 

SO 

*83 

569 

11 

3 

V.  small 

Klunker 

19 

171 

/ 

/ 

V.  small 

FO 

79 

524 

39 

/ 

0 


*:  Error  ^ 0.01  in  off-shock  Cp  and  shock  location 

Rei^.  1 and  3.2  describe  the  extension  of  these  methods  to  lifting  airfoils. 

Here  again  both  FO  and  SO  methods  proved  to  be  of  value.  The  error  in  circulation  of 
the  converged  result  (as  measured  against  a "standard"  result  computed  with  a very 
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large  FD  mesh)  is  plotted  in  Figure  1 as  a function  of  effort  for  each  matching 
method.  It  was  also  discovered  that  the  accuracy  and  cost  of  the  results  was 
indpendent  of  the  location  of  the  circulation  vortex,  an  important  improvement  on 
the  asymptotic  methods. 

Ref.  3.3  describes  the  conclusion  of  the  2D  effort.  Here  the  location  of  the 
matching  boundary  was  moved  to  halfway  between  mesh  points,  and  it  was  proved  that 
the  FO  method  could  converge  to  good  results  even  with  the  lower  "far  field"  boundary 
along  the  lower  surface  of  an  airfoil  (Figure  2),  as  long  as  this  surface  was  sub- 
critical.  The  SO  method  diverged  for  this  case,  and  showed  instability  and  poor 
results  in  other  cases  where  the  matching  boundary  went  through  a subcritical  but  rapidly 
changing  area  of  the  flow  field. 

The  3D  FE  midfield  is  now  operational  and  testbed  results  have  been  generated. 

A rectangular  wing  of  M = 6 with  constant  NACA  0012  section  was  analyzed  by  using 
the  FE  midfield  matching  (Figure  5).  Table  2 gives  the  comparison  between  three  cases; 
large  FD  box,  with  no  matching,  FO  matching  with  panels,  and  double  matching  with  FE's 
and  panels.  In  the  FE  midfield  matching,  the  FO  region  was  reduced  to  a very  small 
box. 


NUMBER  OF 
MATCHING  BOUNDARY 


TECHNOLOGY 

USED 


OUTER  BOUNDARY 
OF  FD  BOX 


FD  MESH  SIZE 
(X  X Y X Z) 


NUMBER  OF 
FD  NODES 


NUMBER  OF 
FE  NODES 


NUMBER  OF 
PANELS 


TABLE  2 


LARGE  BOX 

SMALL  BOX 

0 

1 

FDM 

FDM 

PANEL 

-15  < X IIG 

0 < yTil'f  7 

0 

~'4i.  £ A 1 1 .35^ 

O ^ X 

O ^2 

64  X 28  X 20 

48  X 18  X 14 

35840 

12096 

0 

0 

0 

108 

COMPARISONS  OF  MESH  SIZE 


VERY 

SMALL  BOX 


■'5  1 ;<  1(./V 
^ 3,6 
t .7/ 


42  X 18  X 12 


Note:  Lengths  are  based  on  chord  = 1 
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FE  and  FE  fields  were  iterated  alternately  for  each  far  field  update  and  the 
potentials  at  non-common  points  on  matching  boundaries  were  approximated  by  Lagrangian 
interpolation.  Figure  6 shows  the  results  when  trilinear  elements  were  used  and 
Table  3 demonstrates  a further  reduction  in  computational  time.  The  FE  midfield 
matching  solution,  like  the  first-order  solution,  is  indistinguishable  from  the 
large  box  solution,  while  it  converged  in  about  20%  less  computational  time,  even 
though  some  computational  time  was  used  for  the  Lagrangian  interpolations  between 
FD  and  FE  regions. 

Stability  and  convergence  of  the  method  were  checked  by  testing  several  parameters 
and  the  following  requirements  were  found.  The  linearization  of  the  nonlinear  term  in 
TSDE  should  be  compatible  in  FD  and  FE  algorithms.  The  interface  boundary  should  be 
at  least  one  node  width  outside  of  the  supersonic  FD  area.  Lastly,  in  the  FD  box 
we  used  the  customary  relaxation  factors  of  1.5  for  elliptic  points  and  .9  for 
hyperbolic  points,  but  we  found  that  instability  resulted  from  the  use  of  this 
over-relaxation  in  FE  region  when  alternate  iterations  were  used. 


IV.  CONCLUSIONS 


The  first-order  method  of  far-field  matching,  whereby  the  boundary  conditions 
at  infinity  for  a transonic  flow  problem  with  subsonic  freestream  are  applied  using 
linear  panel  methods  in  the  far  field,  is  a success  both  in  20  and  in  30.  It  allows 
comparatively  small  finite  difference  meshes  to  be  used  without  fear  of  "wind- 
tunnel  -wall  -like"  far  field  errors,  resulting  in  cost  savings  and  increased  confidence. 
(An  estimate  for  the  far  field  error  is  easily  computed,  as  explained  in  ref.  1.). 

The  second  order  method,  successful  in  some  of  the  simpler  2D  cases,  is  touchy 
and  cannot  easily  be  extended  to  30.  Its  alternative,  the  finite  element  mid-field, 
appears  to  be  stable  and  accurate,  and  is  likely  to  yield  cost  advantages  over  the 
first  order  method. 

The  first  order  method  is  not  hard  to  code  and  should  see  immediate  application, 
not  only  in  transonic  small  disturbance  equation  codes,  but  also  in  other  similar 
problems  like  full  potential  equation  codes  and  even  separated-variable  unsteady 
flows  (Helmholtz-like  equations).  The  finite  element  midfield,  being  more  complex 
to  code  and  yielding  less  advantage,  will  probably  still  find  its  way  into  the 
mature  forms  of  these  codes  - a 20%  saving  ought  not  to  be  scorned  in  a high  intensity 
production  use  environment. 


APPENDIX  A 


FINITE  ELEMENT  FORMULATION 


The  Guderley-von  Karman  transonic  small  disturbance  equation,  (1)  is  solved 
in  the  finite  element  midfield  as  well  as  in  the  finite  difference  box  inside. 

In  this  study,  the  finite  element  formulation  is  applied  only  to  an  elliptic 
field  and  an  iterative  process  is  used.  For  a given  Dirichlet  condition,  the 
calculus  of  variations  introduces  a variational  functional 

I (4) 

whose  first  variation  yields  the  TSDE.  That  is,  a function  which  minimizes  the 
functional  is  the  solution  to  the  £f|.(l).  For  non-Dirichlet  (Eq.  4)  type 
boundary  conditions,  a surface  integral  term  must  be  added  to  Eq.  (4).  In  enforcing 
normal  derivative  continuity,  it  proved  to  be  more  convenient  to  update  the  solution 
at  the  boundary  iteratively,  when  an  iterative  process  is  used,  than  to  add  this  term. 

For  the  finite  element  approximation,  a functional  representation  is  used  for 
the  approximate  solution  such  as 

^ (5) 

where  the  are  the  basic  shape  functions,  already  defined,  and  Os  are  the 

coefficients  to  be  estimated.  Since  only  rectangular  elements  are  used  in  this 
study,  a factorization  can  be  used  for  simplicity: 

• Z c.ji.  (6) 

Then,  from  Eq.  (4)  and  the  first  variational  principle 

S;t^^X/atcj  (7) 

the  discretized  finite  element  algorithm  is  derived  by  substitution.  However,  the 
resulting  FE  formulation  becomes  nonlinear  due  to  the  second  term  in  the  RHS  of  Eq.  (4). 
As  in  the  finite  difference  algorithm,  the  nonlinear  term  was  linearized  and  updated 
iteratively.  Hence,  Eq.  (4)  becomes 

xcf)--  ^ 

where  superscript  ^ stands  for  iterate  and  subscript  e is  the  local  element 

in  which  the  velocity  is  evaluated.  It  can  be  shown  that  this  variational  finite 
element  algorithm  in  the  purely  elliptic  field  yields  precisely  the  same  finite 
element  equation  as  the  Galerkin  approach  in  which  tne  sh^e  function  itself  is 
used  for  the  weighting  function. 
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Basic  Shape  Functions  in  the  FE  Formulation 


Only  the  one  dimensional  case  is  shown  in  this  appendix  since  the  3D  shape 
functions  are  factorized  in  the  rectangular  element.  From  Eq.  (5),  the  approximate 
solution  is  defined  by 

(8) 

If  we  define  the  basic  shape  functions  in  each  element  ( J as  shown  in 

Figure  7,  we  have  . ^ 

O)  ' -T.'’  J [ -I  - ' X.  j 

in  a linear  element  and 

for 

in  a quadratic  element. 

Fig.  8 shows  the  difference  in  approximating  a curve  with  these  respective  elements. 

A linear  element  is  an  approximation  of  c/s  which  gives  (<  but 

linear  variation  in  each  element.  Hence  it  is  C®  continuous. 

The  quadratic  element  gives  better  approximation  to  a smooth  function  at  interior 
points  of  the  domain  and  guarantees  C continuity  since  the  derivative  of 
is  continuous  everywhere. 
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L INEAR  KUTTA  CONDITION 
HALF  PARABOLIC  AIRFOIL 
THICKNESS  - .06c 

Moo-  .77 


FIGURE  2 BOUNDARY  MEETS  LOWER  SURFACE  OF  AIRFOIL  IN  MIXED  FLOW 
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